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Studies by video microscopy on fluctuating colloids measure the real-space cross-correlations in particle motion. This set of 
correlations is then treated as a matrix, in order to study the spectrum and mode structure. We show that in general the modes 
are modified by the truncation of the full real-space correlations. We perform a theoretical analysis of the truncation, find the 
boundary conditions imposed by the truncation, and propose practical windowing strategies to eliminate artefacts. We study the 
problem from various perspectives, to compile a survey for experimentalists. 
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1 Introduction 

Many experiments have studied the fluctuations of two and 
three dimensional colloidal solids (both crystalline, amor- 
phous and glassy) using video and confocal microscopy^^]. 
A section of a large sample is observed and recorded; fluctu- 
ations of the subsystem are then analysed off-line. In particu- 
lar one wishes to measures the mode structure and dispersion 
relations of the medium to extract material properties such 
as elastic constants. An important technical question which 
comes up is how to perform analysis on data which is gath- 
ered within a window which is smaller than the full sample. 
One of the main reasons for working with such truncated data 
is of course the hope of eliminating the uncontrolled influence 
of walls by imaging well within the sample using confocal 
techniques. We will show, however, that truncating the full set 
of fluctuations outside of the observation window introduces 
effective boundary conditions which rather unexpectedly lead 
to errors in measurements of material properties. 

The crucial question is how the result of the observed mode 
structure converges when the observation volume increases. It 
is generally assumed that convergence is assured for large sys- 
tem size. After all, Weyl's theorem on the density of states of a 
vibrating system shows that the density of states is asymptoti- 
cally independent of the system shaped thus theorists freely 
interchange the true boundary conditions of a physical sys- 
tem by a mathematically convenient choice of periodic modes. 
However, we show in this paper that the assumption is danger- 
ous when using Fourier analysis of observed amplitudes. To 
be concrete let us now summarize the kinds of data and anal- 
ysis which are available: 

Video or confocal observation of colloidal solids generates 
data sets consisting of the position of particles recorded over 
many frames. One then calculates the mean position of each 
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particle and the correlations in the displacements, w;(r), by 
evaluating the two-particle function 



C i j(r,r') = (u i (r)uj(r')) 



(1) 



with directional indexes i, j, and with r being the reference 
positions of the N particles. The average is over the recorded 
frames. We work in d-dimensional space where the interest- 
ing values are d = 2, 3 when imaging colloids at interfaces, or 
within volumes. 

Three methods of analysis suggest themselves: 

(A) diagonalisation of the correlation matrix eq. ([]]) of dimen- 
sion (dN) x (dN) in order to study its eigenvalues and 
eigenvectors, 

(B) evaluation of amplitudes in Fourier spac^] 

"i(q)=L^ q - r "«(r), (2) 

r 

followed by a similar study of the matrix 

C ij (q,q , ) = (u i (q)uj(q!)), (3) 

(C) direct use of only the diagonal Fourier coefficients 

CiM,-q) = (ui(q)uj(-q)) (4) 

which reduces, for a Bravais lattice, to a set of N matrices 
of size d x d. 

The full diagonalisation of the matrix Qj allows one to plot 
visually seductive pictures of mode structure in the truncated 
system. The question arises as to the exact link between these 
plots and the standard treatment of modes in an elastic medium 

* We differentiate in notation between the reciprocal vectors of the full (possibly 
infinite) system, k, and q for those of the windowed system. 
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via the wave equation. In particular are there effective bound- 
ary conditions introduced by the truncation which influence 
the modes? In this paper we show that in a simple square 
geometry eigenmodes of the correlation matrix do not have a 
pure longitudinal or transverse nature, unlike a bulk sample. 
We quantify this effect in a study of a two-dimensional elastic 
medium. We also note that the decomposition into longitudi- 
nal and transverse is in ambiguous in finite geometries. 

The second formulation in terms of a matrix of Fourier co- 
efficients requires some care as to the definition of the eigen- 
value problem in order to be equivalent to the real space 
formal, in particular in a disordered system plane waves do 
not form an orthonormal basis; thus one must study a gener- 
alized, pencil eigenvalue problem. We will not consider this 
case further in the present paper. 

The use of the third method, studying Qj(q, — q) seems 
particularly practical in large experimental systems because 
it avoids the expensive diagonalisation that is required for the 
other two cases. It might be expected to give exact results for 
elastic moduli in crystals, and good estimates in disordered, 
non-glassy solids. However, we will show numerically that 
Qj(*l> ~ q) i s contaminated by the truncation of the full corre- 
lation functions. We explain the result visually with an anal- 
ogy with the diffraction pattern of the observation window. 
We propose a simple practical solution to the contamination 
by the use of alternative windowing functions with superior 
properties in Fourier space. 

2 Elastic theory 

In this section we resume the results from elasticity theory 
that we will require. In a three-dimensional elastic medium 
the quadratic fluctuations about the energy minimum, as well 
as the propagating modes of a crystal are deduced from the 
Christoffel matrix^. In a cubic solid this has the form 

D ik (k) = X8ij8 k i + n(8i k 8ji + SuSjk) + vS im kjk h (5) 

with Lame constants A, jl and anisotropy v. This expression is 
also valid under uniform, isotropic stress such as the pressure 
which must be applied in non-bound colloidal solids. The free 
energy of small fluctuations about the equilibrium position is 
then given by the functional 



E[u] = £ - Mi -(k)^(-k)A 7 (k). 



(6) 



The Green function of the static elastic problem is then the 
inverse of the Christoffel matrix, 



Dij(k)G jk (k) = S ik . 



(7) 



It describes the response of the medium to static forces, as well 
as correlations in position fluctuations which can be measured 



in microscopy. In a face-centred cubic crystal with nearest- 
neighbour central potentials ji = A = — v, see Eq. (12.7) of 
Ref. [T3] Face-centred hard- sphere systems and real exper- 
iments have non-linearities that slightly modify the relation 
between these three constants^. 

A similar mathematical structure describes fluctuations in a 
two-dimensional hexagonal crystal, with however v = 0, im- 
plying that the long-wavelength mode structure is isotropic 
with just two types of modes- longitudinal and transverse. 
This theory can be used to study the statistics of colloidal crys- 
tals at interfaces 5 . 

For many of the discussions in this paper the most important 
characteristic of the matrix eq. §5§ is the scaling in k 2 . This 
motivates the study of a scalar energy function 



E[u] 



dr: 



X 



(8) 



Such energy functions eliminate the need for detailed tensor 
analysis and allow one to transpose well known theorems in 
potential theory to our study of truncation artefacts. In partic- 
ular for eq. d8l) we find 



G(k) 



1 



(9) 



Then the real-space Green (in d = 3) function is given by the 
Fourier transform of eq. d9]> 



G{r-v>) = Je^G{k)^ 



dk 



i 



47rA|r-r / | 



(10) 



The correlations are translationally invariant. The correspond- 
ing correlations for isotropic elasticity can be found in stan- 
dard references 15 , again the decay in separation varies as 
l/|r — r'|, with additional tensor structure. Expressions with 
cubic anisotropy are treated in the recent literature 16 . The ex- 
perimental correlation matrices are related via equipartition to 
the Green function: 



(ui(k)uj(-k)) =k B TGij(k), 



(11) 



where the correlation is is evaluated for the full system and not 
truncated to a window. 

Experimentalists also study projected correlations, rather 
than the full three-dimensional problem. The experiments thus 
determine a slice of the full correlation matrix. In this case we 
need to determine the effec tive energy function of the sliced 
system. It can be shown^Elthat projection from d to d — 1 di- 
mensions changes the dispersion law in elastic theory from k 2 
to | where is a wave vector in the projected space and 
k is the wave vector in the starting space. 
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3 Corruption of correlations by truncation 

To demonstrate the problem of working with only the diagonal 



elements Q/(q, — q) in truncated data (method (C) in the intro- 



duction) we here present results of molecular dynamics sim- 
ulations performed in two dimensions on a hexagonal crystal 
of hard spheres using event driven methods^. In Fig.[lJ, we 
have analysed the mode structure of the whole, periodic, sys- 
tem. As noted above the long-wavelength mode structure is 
described by two elastic constants, A, jl and there is rotational 
invariance in long- wavelength correlations. In the figure this 
results in there being just two independent intersects for small 
wave vectors when we plot co/k as a function of k. Here CO 2 (k) 
is defined as the eigenvalues of the 2 x 2 matrix G//(k) for a 
given vector k. For the full simulation volume we conclude 
that we are able to extract the effective elastic properties from 
the diagonal values Q/(k, — k). Eq. (pj} therefore holds. 

However, observation of the data in a finite observation win- 
dow leads to considerable modifications in the result, Fig. [I]}. 
We observe a breakdown in rotational invariance, as shown 
by the splitting of previously degenerate modes at small wave 
vectors. In addition the situation does not improve on increas- 
ing the system size. The main point of the first part of the 
present paper is understanding the result of Fig. [TJ) analyti- 
cally and finding analytic and numerical methods which al- 
low one to restore the correct symmetries to the data, in order 
to correctly characterize experimental systems and simulation 
data. 

In Fig. [2] we have performed similar analysis of a three- 
dimensional simulation. The results in panel (a) are more 
complicated than in two dimensions because of the influence 
of cubic anisotropy. But one clearly sees that the very highest 
mode in panel (a) is displaced to lower CO in panel (b), leading 
to important modifications in the effective elastic properties 
that one would deduce from the data. 

We must conclude that truncating data in a larger experi- 
mental system does not lead to a satisfactory method for mea- 
suring elastic properties. We now demonstrate analytically the 
origin of the problem. 

3.1 Windowing theory 

Let us window the data of an experiment or simulation to a 
region w and consider the general two-wavevector transform 

Gw(q,q')= / / e-'' (qr+qV) G(r/)rfrrfr', (12) 

J W J W 

where the integrals are over the observation window. We gen- 
eralize by considering a weighting function W(r) which is 
an arbitrary positive function. We will always normalize this 
function such that 

Jw 2 (r)dr = J\W(k)\ 2 ^ = l. (13) 
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Figure 1 Windowing artefacts in two-dimensional data, (a) Molec- 
ular dynamics data for a hexagonal crystal, analysed for its mode 
structure in periodic boundary conditions. For small wave vectors k 
we see just two branches corresponding to longitudinal (upper, dark) 
and transverse (lower, light) modes, (b) The data is extracted from 
an observation window, half the system size in both directions, anal- 
ysed in the same way, plotted are o)/(q)/|q|. The modes split and 
symmetry-related directions are no-longer equivalent. Note windows 
have non- square shape adapted to the underlying hexagonal lattice 



We find in Fourier space that 

r dk 
Gw(q,q ; ):= / W(q-k)W(k + q')G(k) j-^. 



(14) 



Windowing to a cube corresponds to a constant W(r) within 
the observation zone so that W(k) can be expressed in terms of 
the sine function, well known from the theory of diffraction. 



W(k) 



L d l 2 {\smc fklL 

i=l 



(¥)■ 



(15) 



where the product is over the number of dimensions. This type 
of window will be analysed in detail below. 

We would hope that in the limit of large windows the matrix 
G w (q, — q') is dominated by its diagonal elements and resem- 
bles in some sense G(k), thus W(k) should act as a 5-function. 
We now examine the conditions in which this happens. To un- 
derstand the failures of Fig.fTJ) let us look at an estimate for the 
diagonal elements G w (q, — q). We consider the scalar theory 
withG= 1/Xk 2 . Then 



Gw(q,-q) = / \W(k-q) 



1 dk 



xk 2 (2%y 



(16) 



|W(k)| 2 is simply the Fraunhofer diffraction intensity of the 
measurement window^] 

f The cosine or sine transforms of the correlations can similarly be expressed 
in terms the combinations [G w (q, — q) ± G w (q, q)]. 
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Figure 2 As for Fig. [T] but for three-dimensional analysis. While 
the modification of the mode structure is less drastic than in two di- 
mensions the longitudinal structure is significantly different in the 
original and windowed data sets. 



3.1.1 Gaussian window Explicit progress can be made 
using (idealized) Gaussian windows with 



1 



W(k) = (4a 2 7i 



-r 2 /2a 2 



-a 2 k 2 /2 



(17) 
(18) 



where a is a measure of the real- space width of the window. 
Potentially large contributions to the integral in eq. ( fT6] ) can 
come from the central peak of W(k— q) for cr|k — q| = 
or from the divergence of G(k) at k = 0. However the very 
rapid decay of W(k) implies that the contribution of the in- 
tegral from k = is very small if oq ^> 1. We perform the 
integral in eq. ( p"6| ) with a Gaussian window and find 



G w (q,-q)=e~ 



sjna 
Xq 



erfi(<7g), 



(19) 



With erfi the imaginary error function. When z is small, 
erfi(z) & 2z/\ /f 7l, implying that 



G w (q,-q) : 



for oq small. 



(20) 



When z is large e z erfi(z) ~ 1/ Vn(z 



G w (q,-q) : 



1 



1 



2(oq) 2 



-z~ 3 /2), so that 
for oq large. (21) 



Thus with this well-behaved windowing function G w (q, — q) 
does indeed converge to the desired limit, G(q), where cor- 
rections are higher order in l/aq. While for very small values 
of aq the Gaussian window leads to an overestimate of the 




Figure 3 Fourier transform, |W(k)| 2 , of the weighting function for 
a square sampling window, eq. (T5J. Perpendicular to each side of 
the window the asymptotic decay of the function is slow, l/k 2 . The 
decay in general directions is faster, l/k 4 . The slow decay along 
the symmetry directions leads to major artefacts in the reconstructed 
spectrum. 



elastic modulus, for most reasonable values of the parameter 
we expect a small underestimate in the elastic modulus. The 
theory of windowing with a tensorial Green function is given 
in Appendix |A| where we show that we obtain the correct an- 
swer for the both the longitudinal and transverse modes in an 
isotropic medium when aq large. 

3.1.2 Discontinuous windows We now return to the 
windowing function eq. ( p3] ) which corresponds to truncat- 
ing the data. We plot the diffraction intensity for a two- 
dimensional square in Fig. [3] One sees that the diffraction 
pattern is characterized by a notable "cross"-structure in the 
directions (1,0) and (0, 1). In these two directions the enve- 
lope of Fourier coefficients decays slowly as l/k 2 , while off 
axis the decay is faster, l/k 4 . We now show that the decay in 
l/k 2 leads to convergence of the integral in eq. ( fT6| ) to an in- 
correct value and is the origin of the breakdown in rotational 
invariance in Fig. [T] The splitting that occurs in Fig. [T] in- 
deed corresponds to corruption of modes which are parallel to 
the slowly decaying directions in the diffraction pattern of the 
truncating box. A third mode in the figure, which is equiva- 
lent by symmetry of the hexagonal lattice but not of the win- 
dowing function, remains unaffected by truncation. Note that 
such high symmetry directions are those that are the most nat- 
ural to analyse in an experimental setup. The problem is that 
this diffraction pattern with envelope l/k 2 is not sufficiently 
"close" to a 5-function. Qualitatively we see that we are try- 
ing to extract a signal with a power spectrum in l/k 2 with a 
discontinuous function which has the same power spectrum. 

We now argue quantitatively by estimating the integral 
eq. fT6] ) with the weighting function of Fig. [3] Again two 
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Figure 4 As Fig. [I] but using a Hann window in (b). The windowing 
artefacts are largely eliminated, but the longest wavelength modes 
are modified due to mode mixing, See Appendix A. Skew window 
adapted to simulation box. 




Figure 5 As Fig. [|but using a Hann window in (b). The windowing 
artefacts are largely eliminated, but the longest wavelength modes are 
modified due to mode mixing, compare Appendix |A| Skew window 
adapted to simulation box. 



potentially large contributions come in the neighbourhood of 
k = and k = q. For the integral near k = q, the normal- 
ization of W is designed so that the contribution is exactly 
that corresponding to the physical values of G. It is the ball 
around k = that is particularly problematic, it gives a second, 
non-negligible contribution to the integral. With the natural 
choice of wave vectors qi = 2nn/L the pole at k = is placed 
on a zero of the diffraction pattern, however the contribution 
of a ball of size Ak ~ l/L includes the nearby maxima of 
the diffraction pattern. The envelope of |W(k— q)| 2 is slowly 
varying and can be replaced by a typical value L d / (qL) 2 . In 
the x-direction we have |W(k-q)| 2 - \W (q)\ 2 (k x L) 2 and find 
the contribution the contribution of the ball, 

l/L l/L 

, , m2 f „ , 2 k d ~ l dk L d ~ 2 f k d ~ l dk 1 

o o 

This contribution adds to the &(\/Xq 2 ) contribution near 
k = q. The sine window over-estimates G, and thus naturally 
underestimates effective elastic moduli. In Fig.fTJ) the disper- 
sion curves which are perpendicular to the window limits are 
indeed lower than the correct values. 

3.2 Improved weighting functions 

A detailed study of the spectrum properties of windows in sig- 
nal processing is given by NuttallQ^. We see that the (unre- 
alisable) Gaussian window has excellent properties due to the 
suppression of the contribution of the integral at k = 0. Sinc- 
like windows are inadequate due to their slow decay in Fourier 
space. We now show that a product of Hann windows 

W(x) = (1+cos2^jc/L), -L/2<x<L/2, (23) 



has superior Fourier properties. We find 

W(k)~2 sinc(£L/2) + 

sinc(£L/2 -n) + sinc(/:L/2 + n) , (24) 

for which at large wave vectors |W(&)| 2 ~ l/k 6 . Such a 
rapidly decaying function behaves as a 8 -function when tested 
against a Green function in l/k 2 . 

We reanalysed our two-dimensional data in Fig. [4] using 
this Hann window and find greatly improved results, including 
restoring of rotational invariance in the data. If one repeats the 
analysis of contributions to G w coming from k = one sees 
that they are now sub-dominant. One can trust the results to 
reconstruct the true dispersion law. 

The corresponding results in three dimension are given in 
Fig. [5] We note that results for the longitudinal modes drop 
sharply for small q, while the transverse modes rise. We give a 
theory of this effect in Appendix[A|where we analyse the prob- 
lem of Gaussian windowing in three dimensional isotropic 
media. We show that it is due to a mixing of longitudinal and 
transverse modes which occurs for small wavevectors and cor- 
rupts the dispersion relation. Note this is an effect which only 
occurs for a few low-lying modes- in distinction to the prob- 
lems with sharp windows which lead to all modes in certain 
directions being corrupted. 

3.3 Observation of slices and projection 

What is the situation in projection geometries where the ef- 
fective dispersion law is G(k) ~ l/\k\, cf. Refs. 161171 ? Here 
the use of a sinc-function indeed reproduces the wanted lead- 
ing term in l/\q\ with the correct pref actor, the integral near 
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Figure 6 Projected dispersion curves from the simulation in Fig. [2] 
The projected system is two-dimensional and isotropic, (a) full sys- 
tem, (b) truncation window, (c) window with Hann weighting func- 
tions, see eq. J23J. 



k = q is however less divergent than in the case treated above, 
so that the origin contributes a sub-leading correction in l/q 2 . 
At least asymptotically the correct dispersion relation is ob- 
served in G w (q, — q), though with a leading correction which 
can no-doubt be improved with the use of a windowing func- 
tion which falls to zero at the edge of the observation zone. We 
did indeed simulations to test this point and found the correct 
reconstructed dispersion law, Fig. [6] 

3.4 Off-diagonal elements 

For the ideal case of the Gaussian window one can estimate 
the off-diagonal equivalent to eq. ( [2T] ). Via a saddle point cal- 
culation we find 



G w (q,-q') 



(q + q 



A2' 



-cr 2 (q-q') 2 /4 



(25) 



Higher order corrections can also be calculated. An interest- 
ing, but difficult, problem would be to solve for the eigenvalue 
structure of this effective matrix. 



4 Mode structure of the full correlation matrix 

We now turn to study of the full correlation matrix Q^r,^) 



of a truncated system, method |(A)| in the introduction. An ex- 
perimentalist first measures a set of correlation functions and 
assembles them into such a matrix. It is then natural (and 
easy with tools such as Matlab) to study the eigenvectors of 
this matrix. How are the eigenvectors in a truncated system 
related to those of the underlying physical system which is 
described in terms of elastic properties? We wish to relate 
the eigenvalues and eigenvectors of this matrix in the limit 



of small wave vectors to those from a true continuum theory 
in unbounded space. The particular questions that we study 
include the effective boundary conditions for eigenfunctions 
that arise from the truncation. We will make particular use of 
the scalar analogue of elasticity to simplify the analytic cal- 
culations and to display relevant features of the mathemati- 
cal problem. We find the exact analytic solution to the scalar 
problem in spherical and circular geometries but use numer- 
ical methods in square geometries. We find the quantisation 
conditions by direct study of the duplicating property of an in- 
tegral operator, and we show that simple Neumann and Dirich- 
let boundary conditions do not give the correct mode structure. 

4.1 Scalar elasticity 

We have noted that the problem of scalar elasticity is linked to 
the properties of the Laplacian operator which corresponds to 
the Helmholtz eigenvalue equation: 

[V 2 +£ 2 ]yA = 0, (26) 

The correlation matrix from scalar elasticity is then the form 
C(r,r') = l/47r|r-r'|. 

If we construct a (continuum) matrix from the correlations 
we are thus interested in the following integral equation 



L 



1 



v 4;r|r-r / 



-Y(r) d 3 r = Ai//(r'), 



(27) 



The eigenvalue is A, and we are restricted to a finite vol- 
ume V C R 3 . The mathematical difficulty comes from the 
arbitrary choice for the shape of V. 

Consider equation eq. ( |27| ) with r' within the observation 
volume then we can act on this equation with —V 2 acting on 
the y' coordinate to find 



/ <5(r-r')yA(r) d 3 r = yA(r') = -AVV(r'), 
Jv 



(28) 



which is the Helmholtz equation with eigenvalue A = k~ 2 . 
Thus the matrix of correlation functions has eigenvectors \j/ k 
which are closely related to those of the corresponding differ- 
ential equation, however we will now show that the boundary 
conditions are different. We have the equations 



-V 2 G(r-r') = <5(r-r'), 
(V 2 + £ 2 ) % (r)=0. 



(29) 
(30) 

Multiplying eq. ([29]) by %(r) and eq. ([30]) by G(r - r') leads 
to 

J dr [xi/ k V 2 G-GV 2 xif k -k 2 xi/ k G] = -y/*(r'). (31) 
Using Green's second identity we find 

/ dS r -[\i/ k VG-GV\i/ k ] = -\i/ k (r , )+k 2 [ Jr^(r)G(r,r / ). 

Jdv Jv 
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Figure 7 Amplitudes of modes found in a square observation window 48 x 48 pixels, using scalar elasticity. Top left to bottom right decreasing 
eigenvalues of the integral operator. Note that first mode is not constant within the observation volume despite it being the most uniform mode. 
The second mode displays a noticeable non- sinusoidal form. 



We recognize the right-hand side as the eigenvalue equation in 
integral form, thus the condition 

jf dS r -[vft(r)VG(ry)-G(ry)Vv*(r)]=0 Vr' (32) 

is required for % being an eigenmode both of the elastic prob- 
lem ( [26] ) and of the truncated integral equation ( [27] ). For any 
eigenvector this is a set of integral conditions, true for 
each interior point r'. In spherically symmetric domains V, 



see Sec. |4.3| below, they reduce to boundary conditions valid 
on the surface. In the general case, this cannot be assumed to 
be the case. 



with <9nu the normal stress and djyG the normal derivative of 
the Green function, 



:= QjkidjNi. 



(36) 



The truncated and the elastic eigenvalue problems then have 
the same eigenfunctions if the integral relation 

j> dS[(d N G)-u-G-(d N u)]=0 (37) 

holds for any eigenvector u and for any point r 7 (the integrals 
and derivatives in eqs. ([33])-([37]) act on r). 



4.2 Tensor elasticity 

While performed in the simplest scalar form, the above theory 
can be easily replaced by its tensorial equivalent for an elastic 
medium. The only difference is the use of Betti's identity for 
the stress tensor rather than Green's second theorem^S! Let 

A** := Qjkidjdi in general, and (33) 

A* = (A + ju) graddiv +\l V 2 if isotropic, (34) 

then 

<f dS[(d N G)-u-G-(d N u)] = [ [(A*G)-u-G-(A*-u)]dr 

(35) 



4.3 Scalar spherically symmetric solutions 

We now derive the exact quantization conditions in rotation- 
ally symmetric geometries, again studying the scalar vibra- 
tional problem: Regular solutions to the Helmholtz equation 
in spherically symmetric geometries can be written in terms 
of spherical Harmonics Y™ and spherical Bessel functions 7'/, 



*T(r)=7 Z (tr)ir(0)- 



(38) 



£1 = (0, 0) is a solid angle in spherical polar coordinates. For 
/ = we note that 



li/O 



jo(kr) 



sin(Ar) 
kr 



(39) 



1-(TT]|7 



Let us show that m™ is also a solution to the truncated integral 
equation ( [27] ) for correct choices of k. 

As domain we consider the ball of radius R around the ori- 
gin, V = Br(0), so that the integral operator acting on the trial 
function is given by 



J A7Z\y-y'\ 1 V J 



(40) 



Br(0) 



We use the identity (r< = min(r, r') and r> = max(r, /)) 



r-r 



! /=o r > Z ' "•" 1 m=-/ 



to break the integral into a radial and an angular part, and then 
use the fact that after the angular integrals only a single spher- 
ical harmonic survives so that we must evaluate 

y m rR 1 r l 

/= 27TU \^-\ Mkr) ^ r2dr - (41) 

We explicitly split the integration to find 

Y m / rr* r l + 2 rR r' 1 



Y m ( r r /+2 f R r' 1 \ 

We now use the identities^ 

Jz n+l jn-l(z)=Z n+l jn(z), (42) 

/ ^Jn + i(z) = ~Mz), (43) 

jn+1 (Z) + jn-l(z) = <y2n+ ^ jn(z), (44) 



to transform the integrals into 



l- k2 \JlKkT) 2/+1 {jcR y_ 



i ; ' 



(45) 



This is of the original form if ji_i(kR) = which serves 
to obtain possible values for k: 

1 = 0: cos(kR)=0 => kR= {n+l/2)n, (46) 
1 = 1: sin(kR) =0 kR = nn, etc. (47) 

The condition z ji-i (z) = can be reformulated asEl 

z^-Mz) + (l+l)Mz) = 0, with z = kR, (48) 
az 

which is a mixed (Robin) boundary condition for on the 
sphere. Now we understand how we change the original elas- 
tic problem ([26]) if we truncate the domain: Truncating is 



equivalent to applying the boundary conditions ( [48] ) to the dif- 
ferent solutions of eq. p6| ). Note that for different / we have 
different boundary conditions. 

A very similar expansion in (r</r>) can be performed 
in two dimensions for a disk, where the eigenfunctions 
e im6 J m (kr) give rise to the eigen-equation J m -\(kR) = for 
m > 0. For m = the corresponding equation is 



/o(z)+dog(z/z )/i(z) = 0, (49) 
with z = kR and zo a reference radius. 

4.4 Square geometry 

We were unable to solve the integral equation eq. ( [27] ) in the 
geometry of a square or cubic box. We thus proceed by numer- 
ical investigation. We generate the Green function in a large, 
periodic box of dimension 500 x 500 using fast Fourier trans- 
forms 17 . We then truncate the Green function to a window 
of dimension 48 x 48. This Green function is then diagonal- 
ized using standard dense algebra packages included in Mat- 
lab. The structure of the modes that is found is demonstrated 
in Fig. [7] The simplest mode (top left panel) has a largely 
constant amplitude over the observation region, it does have a 
gentle peak at the centre, however, which is some 30% higher 
than the value of the function in the corners. It is interesting 
to note that this mode has finite energy- it does not correspond 
to a zero mode of the Helmholtz system. The next mode (top 
right panel), which is two-fold degenerate, has the nature of a 
wave within the box. However we see that it is definitely non- 
sinusoidal and the contours of constant amplitude bow out at 
the edge of the observation volume. 

The methods generalize to vector elasticity in two dimen- 
sions. We use a modification of the method of Ref. [T7] to 
generate a discretised version of the matrix D in eq. ([5]), with 
V = 0. In particular we choose a discretised dispersion relation 



D(k) = fi(4-2cosk x -2cosk y ) ® 

)(^~ 2cos& x sinfcjsinfcy^ 
^ ^'ysink x smky 2 — 2cos£y> 



(50) 



and use the fast Fourier transform to generate the correspond- 
ing real-space form G//(r). We truncate this Green function 
to a square and diagonalize. The dispersion relation differs 
from our previous choice 17 and has the advantage of preserv- 
ing more properties of the continuum elastic theory that we 
wish to study. In particular our previous choice leads to a sub- 
dominant contribution to the G xy (Ax, 0) . In our previous work 
it decays as I /(Ax) 2 for large separations Ax. The new form 



eq. pO} gives zero for this quantity. 

In the leftmost column of fig. [5] we plot the vector displace- 
ment fields of three of the lower eigenmodes in an isotropic 
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Figure 8 Three eigenmodes of the truncated system: Left panels show the displacement field, colours indicate absolute value (black=small, 
yellow=large). The three columns on the right show the Helmholtz-Hodge decomposition into longitudinal, transverse and harmonic parts. The 
blue/red colour shows the divergence (second column) and the rotation (third column) of the displacement field (red=negative, blue=positive). 
A = \i. 



elastic medium, windowed to 48 x 48 in a system of dimen- 
sions 500 x 500. Both the colour coding in black/yellow and 
the size of the arrows indicate the absolute value of the dis- 
placements. In these modes we find similar bow-shaped struc- 
tures in the modes to those found in the scalar problem. It is 
interesting to note a number of properties of the figures. When 
we look at the distribution of amplitudes in the top panel, the 
maximum occurs at some distance from the edge of the sam- 
ple and displays a vertical gradient in colour - this shows that 
despite the mode being largely longitudinal in nature it does 
display both longitudinal and transverse characters. 

The longitudinal and transverse nature of these modes is 
better studied in an explicit Helmholtz-Hodge decomposi- 
tion, which is displayed in the second to fourth columns of 
fig. [8] This requires some explanation. We want to split a 
vector field into a rotation-free part, and a divergence-free 



part >i/a, with two scalar fields (f) and y/. The vector differen- 
tial operator > := (— d yi d x ) T replaces the curl in two dimen- 
sion^ The rotation-free part should be "longitudinal", and 
the divergence-free part "transverse". However, this decom- 
position is unique only in infinite space. In the present case 
with a finite window, a third component h may be required, 
having zero divergence and rotation. The decomposition then 
reads 

u = V0+>i/A + h. (51) 

After scalar multiplications with V and with >, one finds that 
the scalars (j) and \j/ satisfy Poisson equations with divergence 
and rotation of the original field as sources, 

A0 = V-u, Ay = >-u. (52) 



$ \\f can be considered the z component of the vector potential. 
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j>| 2 /|V0| 2 J\> ¥ \ 2 /|h|' 



1.00 
1.00 
1.00 



0.42 
0.30 
0.02 



0.04 
0.29 
1.06 



0.25 
0.19 
0.42 



Table 1 L2 norms of the vector fields displayed in fig.[8j in the same 
order. 



We now see that in the finite window the solutions to these 
equations are no longer unique but depend also on the bound- 
ary conditions we apply. One may obtain a possible decompo- 
sition into divergence-free and rotation-free parts by solving 
for (j) with some boundary conditions imposed, e. g. Dirich- 
let or Neumann; the rest, u — is then representable as >l/A. 
If specific boundary conditions are required also for y/, then 
there is a third contribution h which obeys the equation 



> h = V h = 0, Ah = 0. 



(53) 



In our case we lack a reasonable justification of such a bound- 
ary condition for (j) or they give rise to visually unreason- 
able fields- for instance normal fluxes which are forced to zero 
at the edge of the box. We rather choose to come as close as 
possible to what we associate with "longitudinal" and "trans- 
verse" in the infinite system. We therefore do not impose a 
boundary conditions at the box edges but extend eqs. ( |52| ) to 
infinity by padding the right-hand sides with zero and requir- 
ing the solutions to vanish at infinity. The solution to this prob- 
lem has already been discussed in s ec.|4.1| and is the convolu- 
tion of the right-hand sides of eq. (|52[) with minus the scalar 
Green function. The scalar fields (j) and \j/ contain the diver- 
gence and the rotation of the original field u, and we believe 
that these solutions are the least perturbed by any boundary ef- 
fects. The harmonic field h is then simply what remains after 
subtraction of and \>\jf. 

The three contributions of eq. ( [5T] ) are displayed in the three 
rightmost columns of fig. [8] (in the same order). The blue/red 
colour coding shows the right-hand sides of eqs. ( |52| ), red 
for negative values and blue for positive. In the second col- 
umn, displaying V0, we see how the arrows all head to the 
minima/maxima of the divergence field. In the second col- 
umn, displaying >i/a, the arrows turn in positive or negative 
sense around the minima/maxima. The harmonic contribu- 
tion in the third column has zero divergence and zero rotation. 
The arrows in all panels are scaled independently to best vi- 
sualize the fields. Fo r quantitative comparison we give the 
L2 norms in table TJ[ One sees that the mode in the first row 
is mainly longitudinal, but that there is an important harmonic 
(quadrupolar) contribution. We found another mode of simi- 
lar symmetry which was dominated by its transverse part (not 

5 The decomposition is not orthogonal so that the sum of the three individual 
terms is not unity 



shown). The last line of the figure shows again a mainly trans- 
verse mode, but here the harmonic field is not quadrupolar 
but is of lower order than the dominant transverse part. Most 
interestingly, the middle row shows a mode where both lon- 
gitudinal and transverse contributions are equally important. 
This shows clearly that the truncation mixes these two natures 
of the modes. 

A numerical study of the scaling of the mode energies with 
size of the truncating box, t confirmed that, as expected, 
eigenvalues scale as £ 2 ; thus A/£ 2 is the object which contains 
information about material properties. However we note that 
eq. ( [49] ) implies that in a disk geometry there can be a slow 
logarithmic cross-over for certain modes. One might expect 
similar logarithmic corrections in the square geometry too. 

5 Conclusions 

Data analysis, with both numerical and experimental data of- 
ten require windowing. When correlation data is simply trun- 
cated it can lead to substantial artefacts in the measured am- 
plitudes and can mislead as to the exact values of elastic con- 
stants. These errors can be considerably reduced by using 
windowing functions which decay faster in Fourier space. In 
particular we found that the Hann window gives good results. 

The experimental analysis of windowed data also gives rise 
to interesting questions as to the nature of the observed eigen- 
modes. We have shown that the correlation functions give rise 
to problems which satisfy an interesting integral condition in- 
volving the correlation functions of the experimental system. 
We have made a study of the eigenvalue problem for scalar 
elasticity and showed how to find exact analytic solutions of 
the integral equation in spherical and circular geometries. In 
square geometries we exhibited eigenfunctions which deviate 
noticeably from plane waves. We note that the use of integral 
equations to characterize observed experimental correlations 
is known in field such as statistical analysis and atmospheric 
physics^] 

A Elasticity with Gaussian window 

We present here the calculation for a tensorial Green function 
in an isotropic three-dimensional medium, when analysed us- 
ing Gaussian windowing. The longitudinal part of the Green 
function of the original system can be written in the form 



G«(k) 



1 



|k><k| 



A +2^ k 2 



(54) 



with unit column vectors |k) = k/k and their adjoint row vec- 
tors (k|. The Gaussian weighting of eq. (16) involves only a 
single external vector quantity q. Thus we can deduce that 

(A+2 M )G^(q,-q) =A(q)I + B(q)\q)(q\ (55) 
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for two as yet unknown functions A and B. We note that this 
form involves both longitudinal and transverse parts in the new 



variables. We take the scalar product of eq. (54) and eq. (55 ) 
with q and secondly study the trace of these equations to de- 
duce that 



A(q)+B(q): 



l^q-k)!^^ (56) 



k 2 (2tt) 3 

3A(q)+B{q) = J \W(q-k)\ 2 ^^ (57) 

We perform the angular integrals, then recognize the radial 
integrals as being related to imaginary error functions: 



A(q)+B(q) = 



1 



-ay \^rerfi(gg) 



2aq 3 



3A(q) +B(q) = e'^^ erfi(cr^) 
The transverse part of the Green function 



fxk 1 



(58) 
(59) 

(60) 



also gives a contribution which can be expressed in terms of 
the functions A and B: 

LlGtHq-q) = (2A(q)+B(q))I-B\q)(q\. (61) 

The full Green function is then the sum of the contributions of 
eq. ( [55] ) and eq. ( |6T] ). 

When aq is large A « l/(a 2 q 4 ), whereas B{q) « l/q 2 - 
Thus the reconstruction does not mix the longitudinal and 
transverse components of the response which are both found 
correctly. For intermediate values of aq, where A cannot be 
neglected, the transverse and longitudinal modes do mix, in a 
manner similar that we found with the small-wavevector re- 
constructions using the Hann window, Fig. [5] In particular the 
longitudinal stiffness is strongly underestimated. As a specific 
example we take X/ji = 1 and plot in Fig [9] the two effective 
values of ft) as a function of aq. 

Matching the mean squared width of a Hann window to a 
Gaussian gives a 2 = L 2 (l/12— 1/(27T 2 )), giving an approx- 
imate relation between our analytic calculations on Gaussian 
functions and practical windows in the experimental situation. 
For the first mode in a square sample for which ^ = we 
find aq « 1.1. 
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